Climate Indices Notebook

This document reproduces the different steps followed before generating the shiny app.

##Including relevant libraries

library(tidyverse)  # includes ggplot2 and readr

── Attaching core tidyverse packages ───────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(plotly)     # interactive plots

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

##Loading indices file

Mset <- read_delim(file="https://data.mbari.org/products/satellite-derived/global-modes/Mset.txt",comment="%",delim="\t")

Rows: 1363 Columns: 8── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (1): MONTH
dbl (7): YEAR, M1, M2, M3, M4, M5, M6
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Mset$time <- as.Date(paste(Mset$YEAR,Mset$MONTH,rep(15,length(Mset$YEAR))),format="%Y %m %d")
indices <- read_csv(file="https://data.mbari.org/products/satellite-derived/global-modes/climate_indices.csv",comment="%")

Rows: 865 Columns: 11── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (11): year, month, MEI, PDO, NPGO, EMI, Nino3, Nino4, AMO, NAO, ATL3
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

nb_time <- length(indices$year)
indices$time <- as.Date(paste(indices$year,indices$month,rep(15,nb_time)),format="%Y %m %d")

##Plotting MEI & M1

plot(indices$time,indices$MEI,type="l",xlab="Time",ylab="Index")
points(Mset$time,Mset$M1,type="l",col="red")

Testing the use of varname & shading

indexname1 <- "PDO"
indexname2 <- "NPGO"
iok1 <- !is.na(indices[,indexname1])
ipos1 <- indices[,indexname1]>0 & !is.na(indices[,indexname1])
ineg1 <- indices[,indexname1]<0 & !is.na(indices[,indexname1])
x <- indices$time
y0 <- rep(0,nb_time)
plot(x,t(indices[,indexname1]),type="l",xlab="Time",ylab=indexname1)
polygon(c(x[iok1],rev(x[iok1])),c(y0[iok1],rev(t(indices[iok1,indexname1]))),col="skyblue")

points(x,t(indices[,indexname2]),type="l",col="red")

Testing 2 Y-axis

Based on https://www.r-bloggers.com/r-single-plot-with-two-different-y-axes/

indexname1 <- "M1"
indexname2 <- "MEI"
iok1 <- !is.na(Mset[,indexname1])
x <- Mset$time
y0 <- rep(0,length(Mset$YEAR))
par(mar = c(5,5,2,5))
plot(x,t(Mset[,indexname1]),type="l",xlab="Time",ylab=indexname1)
par(new=T)

with(indices,plot(time,MEI,type="l",axes=F,xlab=NA,ylab=NA,col="red",pch=16,cex=1.2))
#plot(indices$time,t(indices[,indexname2]),type="l",axes=F,xlab=NA, ylab=NA,col="red")
axis(side = 4, col="red")

mtext(side = 4, line = 3, indexname2,col="red")
legend("topleft",legend=c(indexname1,indexname2),lty=c(1,1),col=c("black", "red"))

Testing ggplot2 - area plot, double axis, interactive plot

indexname1 <- "M1"
indexname2 <- "MEI"
max1 <- max(Mset[, indexname1], na.rm = TRUE)
max2 <- max(indices[, indexname2], na.rm = TRUE)
Mset[, "varneg"] <- Mset[, indexname1]
Mset[, "varpos"] <- Mset[, indexname1]

Mset$varneg[which(Mset$varneg > 0)] <- rep(0, length(which(Mset$varneg > 0)))
Mset$varpos[which(Mset$varpos < 0)] <- rep(0, length(which(Mset$varpos < 0)))

commonxM <- Mset$time %in% indices$time
commonxI <- indices$time %in% Mset$time
subsetM <- Mset[commonxM, indexname1]
subsetI <- indices[commonxI, indexname2]
iok1 <- !is.na(subsetM) & !is.na(subsetI)
r = cor(subsetM[iok1,], subsetI[iok1, ])
legname <- paste(indexname2, " (r = ", toString(floor(r*100)/100), ")", sep = "")

p <- ggplot() +
  geom_area(data = Mset, aes(x = time, y = varneg), fill = "blue") +
  geom_area(data = Mset, aes(x = time, y = varpos), fill = "red") +
  geom_line(data = Mset, aes(x = time, y = eval(as.name(indexname1))), linewidth = 0) +
  geom_line(data = indices, aes(x = time, y = eval(as.name(indexname2))/(max2/max1), color = legname), linewidth=0.5) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*(max2/max1), name = indexname2)) +
  scale_x_date(limits = c(min(Mset$time)-15, max(Mset$time)+15), expand = c(0, 0)) +
  labs(x = "Time", y = indexname1) +
  theme(legend.position = c(0.15, 0.95),legend.key.width = unit(3,"line"),legend.background = element_rect(fill="transparent"),
        legend.text = element_text(size = 11)) +
  scale_colour_manual("",
                      breaks = c(legname),
                      values = c("black"))
p2 <- plotly_build(p)
#style(p2) %>% layout(legend = list(x = 0.01, y = 1))   # needed because ggplotly removes the legend position
p2